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Abstract— A self-starting A-stable implicit linear multistep block method for approximating stiff initial 
value problems (IVPs) in ordinary differential equations (ODEs) is developed. The construction was carried 
out by pairing a three-step top order method (TOM) and a four-step linear multistep method and shifting 
each equation four times to form a block method that generates approximations on ten grid points 


simultaneously. The implementation of the method on some stiff ODEs confirms its efficiency. 


Keywords— A-Stable, Block Method, Linear Multistep Method (LMM) and Top Order Methods (TOM). 


1.0 INTRODUCTION 
First order ordinary differential equations (ODEs) of the form 


yO=f : 
y(to) = Yo 


represent important mathematical models of real-world phenomena. They are useful tools in the formation 


(1.1) 


of epidemiological models [ 1, 2, 3], dynamic systems [4], chemical reactions, electrical circuits [5] etc. To 
this end, first order ordinary differential equations are applied in different fields such as physical sciences, 


health sciences, social sciences and engineering. 


Over the years, many researchers have concentrated on obtaining solutions to equation (1.1) because an 
nth order ordinary differential equation can be solved by reducing it to a system of first order ODEs which 


are easier to programme. 


some first order ODEs or systems are stiff. “Stiffness occurs in differential equations where two or more 
different time scales of the independent variables on which the dependent variables are changing” [6]. A 


stiff ODE can either be linear or nonlinear. 
In linear problems, stiffness is caused by eigenvalues of large negative values. The degree of stiffness is 


max|a| 
min|A|’ 
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measured using the ratio SR = where A denote an eigenvalue [6]. 
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Although, there are many analytical and semi analytical methods of solving stiff ODEs, especially linear stiff 
ODEs, numerical methods remain valuable methods for dealing with nonlinear problems and for 
investigating and analyzing simple special cases [7]. In recent years, with the invention of software such as 
MAPLE, MATLAB and MATHEMATICA, many researchers have developed A-stable numerical methods for 
handling stiff ODEs. 


Block methods provide very high accurate methods that are absolute stable (A-stable) and circumvent 


Dalghist barrier [8], they were introduced to improve stability of methods [9]. 


Since most block methods are constructed using linear multistep methods, they also provide k — 1 starting 
values and simultaneously generate k approximations [8]. For more on the advantages of block methods, 


the following can be consulted [10, 11, 5, 12, 13]. 


The aim of this paper is to develop a 6th order implicit linear multistep method (LMM) and combine it with 
the top order method (TOM) of order 6 to form a pair which is used to construct a one-step implicit block 


method for approximating stiff ODEs. 


The rest of the paper is organized as follows, Section 2 deals with the formulation and analysis of the 
implicit block method, the implementation is carried out in Section 3, while the discussion and conclusion 


are provided in Section 4 and 5 


2.0 METHODS 

In this Section, we derive a symmetric linear multistep scheme of order (k + 2) = 6, which is combined 
with a top order method (TOM) of order 2k = 6, to form a block method that is self-starting. The pairs are 
shifted forward simultaneously four times to give a set of ten equations which are solved to obtain the 


values of the unknown, Yn, ¥n+15Vn+2) ++» Yn+10° 


2.1 Derivation of a Symmetric Method of Order Six (6) 
We derive a symmetric scheme of order (k + 2) = 6, where k = 4. To this end, consider the polynomial 


function represented by f(x,y, y(X))- 


Let 
n 
r 
Ff (XV) = Pan =) (7) an (2.1) 
k=0 
where 
x= x, + rh,dx = hdr 
and 
X—Xn 
=—— 2.2 
r= (2.2) 


A denotes the forward operator defined by A‘ y, = Ak~1y,4, — A*1yy. 
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Let the interpolating polynomial P, (x,,) be defined by, 


rr 1) ree 2) 
Pu(@tn) = fr + 1A fn + A fy + as, + a 
; 2.3 
rir-1r-2)(r-3 rir—-1)..(r-—n+1 
= DORADA) sae pg p MRD ~( den 
4! n! 
Consider the general well-posed first order initial value problem (IVP), 
dy(%n) 

Fe = Fm VE), YOO) = Yo. (2.4) 


The derivation of the symmetric scheme of order (k + 2), for k = 4 is done in two stages, after which the 


final scheme is obtained by summing the results obtained from the two stages. 
STAGE 1 


Integrating equation (2.4) from x,+43 to X,44 we have, 


I dyn) = | fm Ym) ax (2.5) 


Xn+3 Xn+3 


Using the first five terms of equation (2.3) and substituting in equation (2.5) and further simplification 


gives, 


V(%n+4a) — VXn43) = 








xXn+4 
-1 r-D)(r 2 rr-1)(r-2)(r-3 
i [i+ rah + r(r a apt r( a iy ( dC - )¢ at fa| har (2.6) 
Xn+3 , 
or 
Y(Xn+4) — V(%Xn43) = 
ge ea 4\ 27) 
fn 2 tn 12 tn oy tn 720 tn 5 


Substituting the limits and simplifying the R.H.S. of equation (2.7) gives, 


h {4h + 8h, + Ah + 50h + lh] [3h +5 4h +h t ath Ath, 


1 Juss Juss = B[-Zo ofa t gen fars—gafea tqaghnis tang hes] 28) 
STAGE 2 
Integrating equation (2.4) from x, to X41 gives, 
Xn+1 Xn+1 
[. dy (Xn) = i f (Huy Gea) dx. (2.9) 
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Using the first five terms of equation (2.3) and substituting in equation (2.9) and further simplification 











yields, 
ea r(r-1) (r= - 2-3) 
r(r-1)(r-2)(r- 
YOn41) - yYQn) = = J E + rAf, + in + 4! a*s, hdr 
VOn+1) — YOn) = 
pled nea ee gO Te ig, Oe PO Or 1h (2.10) 
Tf 2 Fra 12 Tn ay ae Tn 720 Tn ‘ 
Substituting the limits and simplifying the R.H.S. of equation (2.12) gives, 
2 1 3 19 4 
hf + Ate pth + 5g a a ha (2.11) 
Simplifying equation (2.11) gives, 
251 11 
GAO 9, = (fu tag haus — gy hui t aephaia aye (2.12) 


Adding equation (2.8) and equation (2.12) results to: 
h 
Yn+4 — Yn+3 Ar Yn+1 — Yn = 50 29 n + 94 fins — 66fn+2 a 94 f n+3 + 29F n+4] (2. 13) 
Equation (2.13) is asymmetric scheme of order k = 4. 


2.2 Six Order Top Order Method (TOM) 
Consider the numerical scheme, 


11 9 9 11 1 9 9 1 
~ 60°” ~ 20> "*1 + 79 nt a 60°73 > hl fn a zo fut ar Zo ln ap Zo fnt3 (2. 14) 


Equation (2.14) is a 3-step TOM formula of order 6 which was considered by [12]. Analysis of the TOMs 


show that they are higher order accurate methods of order p = 2k but are unstable methods. 


Equation (2.14) is combined with the derived symmetric method, equation (2.13) to develop a block 


method for the solution of first order stiff ordinary differential equations. 


2.3 Construction of the Block Method 


The block method is developed by replacing n by n + 3 in equation (2.14) and combining with equation 


(2.13) to form a pair which is shifted four times as presented in block formation. 
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2.4.1 The Order and Consistency of the Block Method 
Order of the Block Method 
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The order of the block method is carried out following the method of block formation and is confirmed to 
be of order 


p=(6, 6. 6, 6, 6, 6 6, 6, 6, 6)? 


with error constant 


271 a7 331 29 317 37 43 37. 2549 ay 
60480 3780 20160 945 12096 1260 1728 945 20160 91476 








Pn+1i = ( 


2.4.2 Convergence Analysis 


Consider the general convergence analysis which is given by the matrix formation of the block method are 


given as 
Ay Your + AoYy = ABoFy + AB, F41 where A, Ap, By and B, are M x M matrices 
Multiplying through by Aj? gives 


Az*AiYn41 + AT*AgYn = RAT*Bo fn + RAT Bi fn4a 
or 


TYn41 + CY¥y = hDfy + RE fn4i (2.15) 
where I = Aj1A,,C = Aj1Ap,D = Aj1By and E = A;1B,. 


Using the test function y = Ay and substituting Ay, for f, and Ayp41 for yn41 leads to 
TYn41 + CY_ — ARE Yn41 — AND yn = 0 (2.16) 


where 


Yn = f Xv Yn), and 


Yns1 = f Ont Ynt1): 


Substitute Z = Ah in equation (2.16) 
TYn41 an CYn = ZEYn+1 = ZD Yn =0 (2.17) 


The characteristic equation associated with (2178) is given by 


(1- ZE)r = (ZD -C) 
or 


r=(1-ZE)1(ZD -C) 

Let r =A 

detUr — A) =0 

where A = (1 — ZE)"1(ZD —C) 


p(r) = detUIr - (1 -ZE)1(ZD - C)] (2.18) 
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Evaluating equation (2.18) gives 


4 AZ BZ C2" DZ EZ? PZ? GL PAZ? IZ jee |. oat 
: KZ10 — LZ9 + MZ8 — NZ7 + OZ® — PZ5 + QZ* — RZ3+SZ2-—TZ+U]} — 09) 
where 
A = 749376414 B = 9072748368 C = 45891413805 


D = 134757262436 E = 281213734995 F = 462030153990 
G = 582582615840 H = 549874580040 I = 367959931200 
J = 516274272000 K = 820643724 L = 9909576702 
M = 50621609595 N = 154860388014 0 = 336857647095 
P = 545142202410 Q = 667471375440 R= 608551131960 
S = 392979211200 T= 161278128000 U = 31755240000 


Recall that Z = Ah, and as h — 0, equation (2.19) reduces to 


U 
9 ee =0 
(r 7) 
or 
r°(r-—1)=0 





OT H=% =3 =% =e =% =] =%y = 7% = Oand yy) = 1 
Since no root has modulus greater than one and |z| = 1is simple, the developed block method is zero stable 


and the block is convergence. 


2.4.3. Region of Absolute Stability 
Equation (2.30) was plotted in MATLAB and the following region confirms that the one-step implicit block 


method is A-stable since the region of stability is the exterior of the circle 
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Figure 1: Region of Absolute Stability of the Constructed Block Method 
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3.0 RESULTS 


In this Section, the block method developed is tested on some first order stiff ordinary differential 


equations. The numerical results, analytical results and absolute errors are displayed in Tables. 


Example 1 
Consider the stiff ordinary differential equations 
y= ly, yO) =f 


Table 1: Numerical/Analytic Results for Example 1 
















































































N Xn Yn y(n) Yn — Y(%n)| Error in [8] 
0 0 1.0000000000 1.0000000000 - - 
1 0.01 | 0.9048374180 0.9048374180 — — 
2 0.02 | 0.8187307522 0.8187307531 9.0000 x 1077° 7.93 x 107? 
3 0.03 | 0.7408182195 0.7408182207 1.2000 x 107? 
4 0.04 | 0.6703200443 0.6703200460 1.7000 x 107? 7.36 x 107? 
5 0.05 | 0.6065306583 0.6065306597 1.4000 x 107? 
6 0.06 | 0.5488116347 0.5488116361 1.4000 x 107? 6.83 x 107? 
7 0.07 | 0.4965853027 0.4965853038 1.1000 x 107? 
8 0.08 | 0.4493289624 0.4493289641 1.7000 x 107? 7.03 x 107? 
9 0.09 | 0.4065696604 0.4065696597 7.0000 x 1077° 
10 | 0.10 | 0.3678794347 0.3678794412 6.5000 x 107? gi2 x 10-8 
Example 2 
Consider the system of two stiff equations 
y' = —20y — 19z, y(0) = 2, and 
z' =-19y — 20z, z(0) = 0, 
Table 2: Numerical/Analytic Results for Example 2 
n Xn Yn (Xn) lyn — Y(Xn)| Zn 2Z(Xn) Zn — Z(Xn)| 
0 0 2.0000000000 | 2.0000000000 = 0 0 = 
1 0.01 | 1.667103289 1.667106708 3.4190 —0.3129963788 | —0.3129929592 3.4196 
x 1076 x 1076 
2 0.02 | 1.438598501 1.438604685 6.1840 —0.5217988451 | —0.5217926620 6.1831 
x 1076 x 1076 
3 0.03 | 1.280806104 1.280812475 6.3710 —0.6600849632 | —0.6600785922 6.3710 
x 1076 x 1076 
4 0.04 | 1.170919276 1.170925510 6.2340 —0.7506596019 | —0.7506533680 6.2339 
x 1076 x 1076 
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5 0.05 | 1.093499889 1.093503496 3.6070 —0.8089589602 | —0.8089553529 3.6073 
x 107° x 107 
6 0.06 | 1.038089291 1.038092172 2.8810 —0.8454397756 | —0.8454368954 2.8802 
x 1076 x 1076 
7 0.07 | 0.9976117888 | 0.9976131096 1.3208 —0.8671758508 | —0.8671745302 1.3206 
x 107° x 1076 
8 0.08 | 0.9672705820 | 0.9672735148 2.9328 —0.8789621102 | —0.8789591780 2.9322 
x 107° x 1076 
9 0.09 | 0.9438328986 | 0.9438280997 4.7989 —0.8840294720 | —0.8840342709 4.7989 
x 1076 x 1076 
10 | 0.10 | 0.9250594058 | 0.9250793294 1.9924 —0.8846154292 | —0.8845955066 1.9923 
x 10-5 x 10-5 
Example 3 


Consider the system of three stiff equations, where t denotes an independent variable while x,y and z 
represent dependent variables 

x’ = —20x — 0.25y — 19.75z, x(0) = 1, 

y' = 20x —20.25y + 0.25z, y(0) = 0, and 

z’ = 20x —19.75y — 0.25z, z(0) = -1, 


Table 3: Numerical /Analytic Results for Example 3 





th Xn x(ty) Yn y(tn) Zn z(t,) 

0 1.0000000000 | 1.0000000000 0 0 —1.0000000000 | —1.0000000000 
0.01 | 0.9800395009 | 0.9800399088 | 0.1776293511 | 0.1776292614 —0.8173831281 | —0.8173832178 
0.02 | 0.9342444613 | 0.9342452012 | 0.3168391390 | 0.3168395534 —0.6732106947 | —0.6732102802 
0.03 | 0.8739732570 | 0.8739740441 | 0.4210195506 | 0.4210202550 —0.5640923891 | —0.5640916846 
0.8077883452 | 0.8077890232 | 0.4947374236 | 0.4947385192 —0.4854612499 | —0.4854601540 
0.05 | 0.7418176353 | 0.7418179491 | 0.5430509762 | 0.5430518386 —0.4322589359 | —0.4322580734 
0.06 | 0.6801550359 | 0.6801551848 | 0.5710142936 | 0.5710151266 —0.3994312400 | —0.3994304070 
0.07 | 0.6252639250 | 0.6252638546 | 0.5833499086 | 0.5833504732 —0.3822555078 | —0.3822549432 
0.08 | 0.5783522634 | 0.5783522936 | 0.5842467243 | 0.5842475754 —0.3765427150 | —0.3765418638 
0.09 | 0.5397093301 | 0.5397085794 | 0.5772654202 | 0.5772648331 —0.3787320615 | —0.3787326487 
































ol a] NI a] wa a] wl wn] B] of] Bg 
o 
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= 















































10 | 0.10 | 0.5089832982 | 0.5089850496 | 0.5653008826 | 0.5653043996 —0.3859285425 | —0.3859250248 
Table 4: Absolute Errors for Example 3 

n th IXn — x(tn)| lYn — V(tn)| Zn — 2(t,)| 

0 0 - al - 

1 0.01 4.0790 x 1077 8.9700 x 1078 8.9700 x 1078 

2 0.02 7.3990 x 1077 4.1440 x 1077 4.1450 x 1077 
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3 0.03 7.8710 x 1077 7.0440 x 1077 7.0450 x 1077 
4 0.04 6.7800 x 1077 1.0956 x 107° 1.0959 x 107° 
5 0.05 3.1380 x 1077 3.6240 x 107° 8.6250 x 1077 
6 0.06 1.4890 x 1077 8.3300 x 1077 8.3300 x 1077 
7 0.07 7.0400 x 1078 5.6460 x 1077 5.6460 x 1077 
8 0.08 3.0200 x 1078 $5110 % 10-7 8.5120 x 10” 
9 0.09 7.5070 x 1077 5.8710 x 1077 5.8720 x 1077 
10 0.10 1.5714 x 107° 3.5170 x 10~° 3.5177 x 10~° 
4.0 DISCUSSION 


A 4th-step linear multistep method (LMM) was developed as presented in equation (2.13). The method is 
symmetric, implicit and zero stable but certainly not A-stable. The solution of its first characteristic 
polynomial in MAPLE reveal that the method has two roots of modulus 1. This shows that the method is 
weakly stable. On the other hand, equation (2.14) is one of top order methods (TOMs) which were 
considered as unstable but highly accurate methods. Ordinarily, combining a TOM and a weakly stable LMM 
cannot produce a method which is capable of approximating stiff ODEs. However, [ 8 and 12] combined the 


TOMs with other LMMs to form block methods which have been successfully implemented on stiff ODEs 


To this end, A pair was formed with equation (2.13), which is a LMM of order 6 and equation (2.14) which 
is a 6th order TOM and used to construct a one-step implicit block method. Analysis was performed on the 


block method and it was verified to be A-stable 


A-stable, self-stating and capable of generating 10 approximations simultaneously. Three examples on a 
single, a system of two and a system of three stiff ODEs were used to test the efficiency of the block method. 
The numerical results showed good approximations of the analytical solutions and also compared 


favourably with other numerical methods. 


5.0 CONCLUSION 
An implicit block method that is self-starting was developed by combining a 4th-step LMM and a 3-step 
TOM. The block method is A-stable. Numerical experiments confirmed that the method provides good 


approximation of stiff ODEs 
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